#include "ecgseekpeaks.h"
#include "ecgfilter.h"

#define WINDOW_SEC (0.090) //(0.160)
#define MIN_RR_SEC (0.220) //(0.3)//(0.4)//(0.200)
#define MIN_RR_SEC_HR (0.200)
#define MAX_RR_SEC (2.0)

#define MAX_VOLT (4.0) //最大伏值
#define MIN_VOLT (0.25) //最小伏值
#define MIN_QS_MS (40.0) //最小QS波宽度
EcgSeekPeaks::EcgSeekPeaks(QObject* parent)
    : QObject { parent }
{
}

void EcgSeekPeaks::recognition(QVector<float>& ecgSignalData)
{
    qDebug() << "进入识别心搏函数";
    int sampleRate = ecg.sampleRate();
    const int MIN_RR = (int)(MIN_RR_SEC_HR * sampleRate);
    // const int MAX_RR = (int)(MAX_RR_SEC * sampleRate);
    //滤基线
    avgFilter(ecgSignalData);

    //低通高通滤波50hz
    // EcgFilter ecgFilter;
    // ecgFilter.initBp(30.0, sampleRate, 10.0);
    // for (int i = 0; i < ecgSignalData.size(); ++i) {
    //     float fdatum = ecgFilter.bp(ecgSignalData[i]);
    //     // fdatum = ecgFilter.hp(fdatum);
    //     // fdatum = ecgFilter.bp(fdatum);
    //     ecgSignalData[i] = fdatum;
    // }
    qDebug() << "滤波完成";
    //寻找R波
    QVector<int> result;
    rPeakPositions(ecgSignalData, result, MIN_RR);
    qDebug() << "寻找r波完成";

    // QHash<int, QHash<unsigned char, int16_t>> beats;
    // int n = 0;
    // for (int i = 0; i < result.size(); ++i) {
    //     beats.insert(result[i], { { 0, 1 } });
    //     n++;
    // }
    qDebug() << "共检出心搏" << beatsModel.data.size() << "个";
    // beatsModel.data = beats;
}
/**
 * @brief 滤基线
 * @param ecgSignalData
 */
void EcgSeekPeaks::avgFilter(QVector<float>& ecgSignalData)
{
    float avgAvlue;
    double sum = 0;
    for (int i = 0; i < ecgSignalData.size(); ++i) {
        sum += ecgSignalData[i];
    }
    avgAvlue = sum / ecgSignalData.size();
    for (int i = 0; i < ecgSignalData.size(); i++) {
        ecgSignalData[i] = ecgSignalData[i] - avgAvlue;
    }
}
//排序后取出百分比所在的值
float EcgSeekPeaks::percentiles(const QVector<float>& data, int percent)
{
    QVector<float> sortedList = data;
    // 对列表进行升序排序
    std::sort(sortedList.begin(), sortedList.end());
    return sortedList[(int)((float)sortedList.size() * (float)percent / 100 + 0.5)];
}

/**
 * @brief 识别r波位置
 * @param signal 心电信号
 * @param result r波位置
 * @param min_rr_width 最小rr间期宽度
 */
void EcgSeekPeaks::rPeakPositions(const QVector<float>& signal, QVector<int>& result, int min_rr_width)
{
    QHash<int, int> allRpeakValues;
    int sampleRate = ecg.sampleRate();
    QVector<int> peakValleyPositions;
    QMap<int, QHash<QString, float>> allRPeaks;
    //检出所有波峰
    findPeaks(signal, peakValleyPositions);

    int currentFragmentStart = -1, highVoltPeak;
    // QList<int> currentMax(10);
    QHash<QString, float> currentMax;

    for (int i = 1; i < peakValleyPositions.size(); i += 2) {
        float value = signal[peakValleyPositions[i]]; //波峰值
        float preValue = signal[peakValleyPositions[i - 1]]; //前一个波谷的值
        float nextValue = signal[peakValleyPositions[i + 1]]; //后一个波谷的值
        float volt = value - std::min(preValue, nextValue); //伏值
        float qsWidth = peakValleyPositions[i + 1] - peakValleyPositions[i - 1]; //宽度
        int preMidpoint = peakValleyPositions[i] - (peakValleyPositions[i] - peakValleyPositions[i - 1]) / 2; //波峰和前一个波谷的中间点
        int nextMidpoint = peakValleyPositions[i] + (peakValleyPositions[i + 1] - peakValleyPositions[i]) / 2; //波峰和后一个波谷的中间点

        QPointF O(0.4 * peakValleyPositions[i] / sampleRate, value);
        QPointF A(0.4 * preMidpoint / sampleRate, signal[preMidpoint]);
        QPointF B(0.4 * nextMidpoint / sampleRate, signal[nextMidpoint]);
        QHash<QString, double> dot = angleDotProduct(O, A, B); //波峰形成的角度
        double angle = dot["angle"];
        double dotProduct = dot["dotProduct"];

        //筛选出符合R波特征的波峰
        if (
            volt <= MAX_VOLT
            && volt >= MIN_VOLT //检测波峰波谷差小于4且大于0.2
            // && qsWidth > (sampleRate * (MIN_QS_MS / 1000.0f)) //宽度大于10ms的波峰
            && angle < 10 //角度小于10
            && volt / (0.4 * qsWidth / sampleRate) > 5 //高宽比大于5
        ) {
            //选出min_rr_width范围内最大波峰作为R波备选
            if (currentFragmentStart == -1) {
                //片段开始
                currentFragmentStart = peakValleyPositions[i];
            } else if (peakValleyPositions[i] - currentFragmentStart >= min_rr_width) {
                if (!highVoltPeak || (highVoltPeak && std::abs(currentMax["rPeak"] - highVoltPeak) > 0.5 * sampleRate)) {
                    //超过最小rr间期，取区间内最大值作为当前R波
                    allRPeaks.insert(currentMax["rPeak"], currentMax);
                    allRpeakValues[currentMax["volt"] * 10]++; //保留一位小数有效值
                }
                //新片段开始
                currentFragmentStart = peakValleyPositions[i];
                currentMax.clear();
            }
            //更新片段内的最大值
            if (currentMax.isEmpty() || value > signal[currentMax["rPeak"]]) {
                currentMax["rPeak"] = peakValleyPositions[i];
                currentMax["q"] = peakValleyPositions[i] - peakValleyPositions[i - 1];
                currentMax["s"] = peakValleyPositions[i + 1] - peakValleyPositions[i];
                currentMax["angle"] = angle;
                currentMax["volt"] = volt;
                currentMax["qsWidth"] = qsWidth;
                currentMax["dotProduct"] = dotProduct;
                currentMax["value"] = value;
            }

        } else if (volt > MAX_VOLT) {
            highVoltPeak = peakValleyPositions[i];
        }
    }
    //计算出现最多次数的伏值
    float normalRpeakValue = 0.0; //推测的R波正常值
    int maxCount = 0;
    for (auto it = allRpeakValues.constBegin(); it != allRpeakValues.constEnd(); ++it) {
        if (it.value() > maxCount) {
            maxCount = it.value();
            normalRpeakValue = it.key() / 10.0f;
        }
    }
    qDebug() << "normalRpeakValue：" << normalRpeakValue;
    //对筛选出的R波做二次验证
    int lastNnIntervals = 0; //最后一次nn间期
    int rrIntervals = 0; //当前跳的rr间期
    int preLabel; //前一跳的标签
    QHash<QString, float> normalNTemp; //标准心搏模板
    QList<int> preRPeaks; //有效心搏
    for (auto it = allRPeaks.constBegin(); it != allRPeaks.constEnd(); ++it) {
        // beatsModel.data.insert(it.key(), it.value());
        // continue;
        // if(beatsModel.data.size() < 30){
        //      qDebug() << "lastNnIntervals" << lastNnIntervals<< "rrIntervals" << rrIntervals;
        // }
        QHash<QString, float> peakValue = it.value();
        QHash<unsigned char, int16_t> beatValue = { { beatsModel.labelIndex, 1 }, { beatsModel.qIndex, peakValue["q"] }, { beatsModel.sIndex, peakValue["s"] }, { beatsModel.test1Index, peakValue["angle"] * 10 }, { beatsModel.test2Index, peakValue["volt"] * 1000 }, { beatsModel.test3Index, peakValue["qsWidth"] } };
        int peakPosition = it.key();
        float volt = peakValue["volt"];
        //第一个心搏的伏值必须在normalRpeakValue的正负0.2之间
        if (preRPeaks.size() == 0 && std::abs(volt - normalRpeakValue) <= 0.2) {
            beatsModel.data.insert(peakPosition, beatValue);
            preLabel = beatsModel.N;
            preRPeaks.prepend(peakPosition);
            continue;
        } else if (preRPeaks.size() != 0
            && std::abs(volt - normalRpeakValue) <= normalRpeakValue * 1) {
            int preRrIntervals = rrIntervals; //前一跳的rr间期
            //当前跳的rr间期
            rrIntervals = peakPosition - preRPeaks[0];
            allRPeaks[peakPosition].insert("rrIntervals", rrIntervals);
            beatValue.insert(beatsModel.rrIntervalsIndex, rrIntervals);
            if (preRrIntervals == 0) {
                //前两跳默认为N
                beatsModel.data.insert(peakPosition, beatValue);
                preLabel = beatsModel.N;
                preRPeaks.prepend(peakPosition);
                lastNnIntervals = rrIntervals;
                continue;
            }
            //判断室早后一跳是否为误报
            if (preLabel == beatsModel.V
                && std::abs(volt - normalRpeakValue) > normalRpeakValue * 0.15
                && rrIntervals < 0.4 * sampleRate
                && (std::abs(peakValue["angle"] - allRPeaks[preRPeaks[0]]["angle"]) > 0.2 * peakValue["angle"]
                    || std::abs(peakValue["volt"] - allRPeaks[preRPeaks[0]]["volt"]) > 0.2 * peakValue["volt"]
                    || std::abs(peakValue["qsWidth"] - allRPeaks[preRPeaks[0]]["qsWidth"]) > 0.2 * peakValue["qsWidth"])) {
                //如果前一跳是V，后一跳形态与它差异大且不符合N的特征，则可认为是V后面的T波被误识别
                continue;
            }
            beatsModel.data.insert(peakPosition, beatValue);
            preLabel = beatsModel.N;
            //rr间期比前一跳提前20%以上的为早搏
            if (rrIntervals < lastNnIntervals * (1 - 0.2)) {
                if (peakValue["qsWidth"] > allRPeaks[preRPeaks[0]]["qsWidth"] * (1 + 0.2)
                    || std::abs(peakValue["volt"] - allRPeaks[preRPeaks[0]]["volt"]) > 0.2 * peakValue["volt"]) {
                    //qrs波宽大20%以上或伏值差距20%以上，为室早
                    beatsModel.data[peakPosition][beatsModel.labelIndex] = beatsModel.V;
                    preLabel = beatsModel.V;
                } else {
                    //房早
                    beatsModel.data[peakPosition][beatsModel.labelIndex] = beatsModel.S;
                    preLabel = beatsModel.S;
                }
            }
            preRPeaks.prepend(peakPosition);
            //计算正常心搏模板
            if (preRPeaks.size() > 2) {
                //生成标准心搏模板
                setNormalNTemp(preRPeaks, beatsModel.data, allRPeaks, signal, normalNTemp);
                if (!normalNTemp.isEmpty()) {
                    lastNnIntervals = normalNTemp["rrIntervals"];
                }
            }
        }
    }
}
/**
 * @brief 寻找所有波峰
 * @param signal 心电信号
 * @param peakValleyPositions 波峰和波谷的位置vpvpvpvpv
 */
void EcgSeekPeaks::findPeaks(const QVector<float>& signal, QVector<int>& peakValleyPositions)
{
    QVector<int> diff_v(signal.size() - 1, 0);
    // 计算V的一阶差分和符号函数trend
    for (QVector<int>::size_type i = 0; i != diff_v.size(); i++) {
        //降低精度
        if (static_cast<int>(signal[i + 1] * 100) - static_cast<int>(signal[i] * 100) > 0)
            diff_v[i] = 1;
        else if (static_cast<int>(signal[i + 1] * 100) - static_cast<int>(signal[i] * 100) < 0)
            diff_v[i] = -1;
        else
            diff_v[i] = 0;
    }
    // 对Trend作了一个遍历
    for (int i = diff_v.size() - 1; i >= 0; i--) {
        if (diff_v[i] == 0 && i == diff_v.size() - 1) {
            diff_v[i] = 1;
        } else if (diff_v[i] == 0) {
            if (diff_v[i + 1] >= 0)
                diff_v[i] = 1;
            else
                diff_v[i] = -1;
        }
    }
    bool firstV = true;
    for (QVector<int>::size_type i = 0; i != diff_v.size() - 1; i++) {
        //从第一个波谷开始存
        int diff = diff_v[i + 1] - diff_v[i];
        if (firstV == true && diff == 2) {
            peakValleyPositions.append(i + 1);
            firstV = false;
        } else if (firstV == false && (diff == 2 || diff == -2)) {
            peakValleyPositions.append(i + 1);
        }
    }
    if (peakValleyPositions.size() % 2 == 0) {
        // 如果是双数，删除最后一个元素
        peakValleyPositions.removeLast();
    }
    // qDebug() << peakValleyPositions.size();
}
/**
 * @brief 计算3个点形成的夹角和点积
 * @param o 中心点
 * @param first
 * @param second
 * @return
 */
QHash<QString, double> EcgSeekPeaks::angleDotProduct(const QPointF& o, const QPointF& first, const QPointF& second)
{

    double dsx, dsy, dex, dey, angle, dotProduct;
    QHash<QString, double> result;

    dsx = first.x() - o.x();
    dsy = first.y() - o.y();

    dex = second.x() - o.x();
    dey = second.y() - o.y();

    dotProduct = dsx * dex + dsy * dey;

    double c = qSqrt(dsx * dsx + dsy * dsy) * qSqrt(dex * dex + dey * dey);
    if (0 != c) {
        angle = qAcos(dotProduct / c);
        angle = qRadiansToDegrees(angle);
    }
    result.insert("angle", angle);
    result.insert("dotProduct", dotProduct);
    return result;
}
/**
 * @brief 根据已识别的心搏设置标准心搏模板
 * @param preRPeaks
 * @param allRPeaks
 * @param signal
 * @param normalNTemp
 */
void EcgSeekPeaks::setNormalNTemp(const QList<int>& preRPeaks, const QHash<int, QHash<unsigned char, int16_t>>& beatsData, const QMap<int, QHash<QString, float>>& allRPeaks, const QVector<float>& signal, QHash<QString, float>& normalNTemp)
{
    if (preRPeaks.size() > 2) {
        // if (normalNTemp.isEmpty()) {
        //     normalNTemp["value"] = 0; //值
        //     normalNTemp["volt"] = 0; //伏值
        //     normalNTemp["qsWidth"] = 0; //qr宽度;
        //     normalNTemp["qsHeight"] = 0; //qr高度差，干扰的情况下qr会有高度差;
        //     normalNTemp["dotProduct"] = 0; //点积
        //     normalNTemp["angle"] = 0; //角度
        //     normalNTemp["rrIntervals"] = 0; //rr间期
        // }
        QList<float> valueList;
        QList<float> voltList;
        QList<float> qsWidthList;
        QList<float> qsHeightList;
        QList<float> dotProductList;
        QList<float> angleList;
        QList<float> rrIntervalsList;
        int rpeak;
        int preRpeak;
        for (int i = 1; i < preRPeaks.size(); ++i) {
            if (valueList.size() >= 10) {
                //取前10跳作为正常模板的基础数据
                break;
            }
            rpeak = preRPeaks[i];
            preRpeak = preRPeaks[i - 1];
            if (beatsData[rpeak].value(beatsModel.labelIndex) == beatsModel.N && beatsData[preRpeak].value(beatsModel.labelIndex) == beatsModel.N) {
                //本跳和前一跳都是N的计入正常模板
                valueList.append(allRPeaks[rpeak].value("value"));
                voltList.append(allRPeaks[rpeak].value("volt"));
                qsWidthList.append(allRPeaks[rpeak].value("qsWidth"));
                float qsHeight = signal[rpeak - allRPeaks[rpeak].value("q")] - signal[rpeak + allRPeaks[rpeak].value("s")];
                qsHeightList.append(qsHeight);
                dotProductList.append(allRPeaks[rpeak].value("dotProduct"));
                angleList.append(allRPeaks[rpeak].value("angle"));
                rrIntervalsList.append(allRPeaks[rpeak].value("rrIntervals"));
            }
        }
        if (valueList.size() > 0) {
            normalNTemp["value"] = averageMedian(valueList);
            normalNTemp["volt"] = averageMedian(voltList);
            normalNTemp["qsWidth"] = averageMedian(qsWidthList);
            normalNTemp["qsHeight"] = averageMedian(qsHeightList);
            normalNTemp["dotProduct"] = averageMedian(dotProductList);
            normalNTemp["angle"] = averageMedian(angleList);
            normalNTemp["rrIntervals"] = averageMedian(rrIntervalsList);
        }
    }
}
/**
 * @brief 计算一个list的中间数据的均值
 * @param listData
 * @return
 */
float EcgSeekPeaks::averageMedian(const QList<float> listData)
{
    QList<float> sortedList = listData;
    std::sort(sortedList.begin(), sortedList.end());
    int removeCount = 0;
    if (sortedList.size() > 4) {
        removeCount = 2;
    } else if (sortedList.size() > 2) {
        removeCount = 1;
    }
    //去掉最大最小值
    for (int i = 0; i < removeCount; ++i) {
        sortedList.removeFirst();
        sortedList.removeLast();
    }
    //计算均值并返回
    double sum = 0.0;
    for (double value : sortedList) {
        sum += value;
    }
    return sum / sortedList.size();
}
